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Abstract. We study the phase behavior of hard spheres confined between two 

parallel hard plates using extensive computer simulations. We determine the 

full equilibrium phase diagram for arbitrary densities and plate separations from 

^4H ' one to five hard-sphere diameters using free energy calculations. We find a 

w , first-order fluid-solid transition, which corresponds to either capillary freezing 

^ , or melting depending on the plate separation. The coexisting solid phase 

consists of crystalline layers with either triangular (A) or square (D) symmetry. 
Increasing the plate separation, we find a sequence of crystal structures from 
• ■ -nA — > (n -|- l)n — > (n -I- 1)A ■ • ■, where n is the number of crystal layers, in 
agreement with experiments on colloids. At high densities, the transition between 
'"^ , square to triangular phases are intervened by intermediate structures, e.g., prism, 

^ ■ buckled, and rhombic phases. 
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1. Introduction 

The physics of confined systems is important in different fields of modern technology, 
like lubrication, adhesion and nanotechnology. The study of simple models is 
instrumental in understanding the behavior of complex systems. As such the hard- 
sphere system plays an important role in statistical physics; it serves as a reference 
system for determining the structure and phase behavior of complex fiuids, both 
in theory and simulations. The bulk phase behavior of hard spheres is now well 
understood. At sufficiently high densities, the spheres can maximize their entropy 
by forming an ordered crystal phase E] . The insertion of a hard wall in such a 
fluid decreases the number of hard-sphere configurations. The system can increase its 
entropy by the spontaneous formation of crystalline layers with triangular symmetry, 
the (111) plane, at the wall, while the bulk is still a fluid [S|. This effect is known 
as prefreezing, and is analogous to complete wetting by fluids at solid substrates. It 
is induced by the presence of a single wall and should not be confused with capillary 
freezing. Capillary freezing denotes the phenomenon of confinement induced freezing 
of the whole fluid in the pore at thermodynamic state points where the bulk is still 
a fluid. This transition depends strongly on the plate separation. The opposite 
phenomenon, called capillary melting, can also occur. The capillary induces melting 
for thermodynamic state points that correspond to a crystal in the bulk. Confinement 
can also change dramatically the equilibrium crystal structure. In 1983 Pieranski 0] 
reported a sequence of layered solid structures with triangular and square symmetry 
for colloidal hard spheres confined in a wedge. The sequence of high density structures 
is determined more accurately in recent experiments ISII^, reporting the observation of 
prism phases with both square and triangular symmetry. Recently Cohen 7 studied 
configurations of confined hard spheres under shear, demonstrating the importance 
of the equilibrium configurations in the rheological properties. Despite the great 
number of theoretical and simulation studies on confined hard spheres 8, 9, 10 , the full 
equilibrium phase behavior is yet unknown. In fact, many of the previous studies were 
based on an order parameter analysis, which fails dramatically in discriminating the 
different structures at high densities and large plate separations. More importantly, 
free energy calculations of confined hard spheres are prohibited so far due to the 
lack of an efficient thermodynamic integration path which relates the free energy of 
interest to that of a reference system, while a further complication arises from the 
enormous number of possible solid phases that has to be considered. Hence, it is 
unresolved whether the experimentally observed phases are stabilized kinetically or 
are thermodynamically stable. 

In this letter, we present a novel efficient thermodynamic integration path that 
enable us to calculate the free energy of densely packed and confined hard spheres, 
with high accuracy close to the fiuid-solid transition. This method allow us to 
determine for the first time the stability of the structures found in experiments. 
To this end we perform explicit free energy calculations to map out the full phase 
diagram for plate separations from 1 to 5 hard-sphere diameters. We report a dazzling 
number of thermodynamically stable crystal structures (26!) including triangular, 
square, buckling, rhombic, and prism phases, and a cascade of corresponding solid- 
solid transformations. In addition, the free energy calculations allow us to determine 
the chemical potential at coexistence, that was unaccessible in previous simulations. 
From the analysis of the chemical potential, we find an intriguing sequence of capillary 
freezing and melting transitions coupled to a structural phase transition of the confined 
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Figure 1. (Color online) Illustration of hard spheres with diameter a, confined 
between parallel hard plates of area A = LxLy and a separation distance H. 



crystal. We note that our new method and results are also relevant for confined 
simple fluids |12lllll IT^ and self- assembled biological systems ^J. In addition, the 
structure of dense packings of spheres explains the shape of for instance snowflakes, 
bee honeycombs, and foams, and it is of great importance for fundamental research, 
e.g., solid state physics and crystallography, and for applications like communication 
science or powder technology |15| . 

2. Model and Method 



Our model system consists of N hard spheres with diameter a, confined between two 
parallel hard plates of area A — L^Ly (Fig. ^|. In each layer we used approximately 
200 particles. We use the packing fraction 77 = 7ra^N/{6AH) as a dimensionless 
density, where H is the distance between two plates. We determine the equilibrium 
phase diagram by performing Monte Carlo (MC) simulations in a box, which is allowed 
to change its shape to accommodate different types of crystals; the ratio L^/Ly may 
vary while H and A are fixed. Trial solid structures are obtained from crystals with 
triangular or square symmetry relaxed with MC moves while slowly increasing the 
density by expanding the spheres. The free energy F for the resulting equilibrated 
structures is calculated as a function of rj and H . We use the standard thermodynamic 
integration technique |16l ITT) , but with a new and efficient path based on penetrable 
potentials, that enable us to change gradually from a non-interacting system to the 
confined hard-core system of interest. The sphere-sphere potential reads 



and the wall-fluid potential 

Vwiiz) = 



eexp(- 
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) if Rij < Uc 
otherwise 

if Zj < crc/2 

otherwise 
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where Rij is the distance between spheres i and j, Zi is the distance of sphere i to 
the nearest wall, A and B are adjustable parameters that are kept flxed during the 
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Figure 2. (Color online) The equilibrium phase diagram of hard spheres with 
diameter cr confined between two parallel hard walls in the plate separation H - 
packing fraction r] representation. The white, yellow and dotted regions indicate 
the stable one-phase region, the two-phase coexistence region, and the forbidden 
region, respectively. 



Table 1. List of intermediate structures In as found in our simulations and in 
the experiments of Fontecha 1^1. 



Phase 


Transition 


Simulation! 


Experiment 


11 


lA - 


-+20 


2B 


2B 


12 


20- 


-»2A 


27^ 


2TI 


13 


2A- 


-*3n 


2Pa + 2Pn 


2Va + 2Va 


14 


SD- 


-+3A 


3n + 3Va +(3B) 


Sn + 3Pn + 3Pa 


15 


3A- 


-^40 


3Pa + 4B 


3-Pn 


16 


40- 


-»4A 


4-Pq -f 47^ + 4Pa 


4-Pn -f 4Pa 


17 


4A- 


-*5n 


4-Pa 


4PA,4-Pn,//§ 


18 


SD- 


->5A 


5-Pn-f4-PA+(5PA)+5-R. 


5-p|| 


19 


5A- 


-+60 


5Pa 


no data 


no 


en- 


■+6A 


5Vn + 5Va 


no data 



simulations, and e is the integration parameter. The hmit e ^ cxi yields the hard-core 
interaction, but convergence of the thermodynamic integration is already obtained for 
e ~ 70k bT. The reference states ( e — Ofc^r) are the ideal gas and the Einstein crystal 
for the fluid and solid phase, respectively. We use a 21-point Gaussian quadrature for 
the numerical integrations and the ensemble averages are calculated from runs with 
40000 MC cycles (attempts to displace each particle once), after first equilibrating 
the system during 20000 MC cycles. We determine phase coexistence by equating the 
grand potentials il = F — fj,N |18| . 
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3. Results 

To validate this approach, we perform simulations of a bulk system of hard spheres and 
we find that the packing fractions of the coexisting fluid and face-centered-cubic (fee) 
solid phase are given by 77/ = 0.4915 ± 0.0005 and rjs = 0.5428 ± 0.0005, respectively. 
The pressure and the chemical potential at coexistence are (3Pa^ = 11.57 ± 0.10 
and /3/i — 16.08 ± 0.10. These results are in good agreement with earlier results 
0E1- Furthermore, to validate the approach for confined systems, we determine 
at bulk coexistence the wall-fiuid interfacial tension /37wfO'^ = 1.990 ± 0.007, and 
the wall-solid interfacial tension for the (111) and (100) planes of the fee phase, 
p-fll^a^ = 1.457 ± 0.018 and f3"/i°°(T'^ = 2.106 ± 0.021. Our resuhs are in agreement 
with previous simulations "20", but the statistical error is one order of magnitude 
smaller due to our new thermodynamic integration path. 

Employing this approach we determine the phase behavior of confined hard 
spheres for plate separations 1 < H/a < 5. Fig. El displays the full phase diagram 
based on free energy calculations in the H ^ rj representation. The white regions of 
the phase diagram denote the stable one-phase regions. The (yellow) shaded regions 
indicate coexistence between fluid and solid or two solid phases, and the dotted region 
is forbidden as it exceeds the maximum packing fraction of confined hard spheres. At 
low densities, we observe a stable fluid phase followed by a fluid-solid transition upon 
increasing the density. The oscillations in the freezing and melting lines reflect the 
(in)commensurability of the crystal structures with the available space between the 
walls. For the crystal phases, we follow the convention introduced by Pieranski ^, 
where nA denotes a stack of n triangular layers, and nD a stack of n square layers. 
For H/a — > 1, the stable crystal phase consists of a single triangular layer lA, which 
packs more efficiently than the square layer. As the gap between the plates increases, 
crystal slabs with triangular (Fig.j^^a)) and square packings (Fig.|2Ib)) are alternately 
stable. We find the characteristic sequence • • • nA — > (n -I- l)n —^ {n + 1)A, which 
consists of an nA — > (n + l)n transformation where both the number of layers and 
the symmetry change followed by an (n -f l)n -^ (n -I- 1) A transformation where only 
the symmetry changes. This sequence is driven by a competition of a smaller height 
of n square layers compared to n triangular layers and a more efficient packing of 
triangular layers w.r.t. square layers. When the available gap is larger than required 
for the nA structure, but smaller than for (n -I- 1)0, intermediate structures may 
become stable. Similar arguments can be used for the intervention of intermediate 
structures in the (n-l-l)n — > (n-l-l)A transformation. Especially at high packing 
fractions, the spheres can increase their packing by adopting interpolating structures. 
In Fig. [21 we report the boundaries of the interpolating regions In. Each region 
represents one or more interpolating structures, that are listed in Tab. ^ according 
to the standard notation. Within the resolution of our simulations, it is difficult to 
draw the phase boundaries of all the intermediate structures in /n, but in Tab. ^ the 
thermodynamically stable structures are listed in the order they appear upon increasing 
H and rj. We also compare our sequence of structures with the experimental one 
jHl. The experiments considered charged particles, but we do not expect that the 
soft repulsion has a strong effect on the observed structures at high densities. The 
agreement is excellent at small plate separations. The buckling phase 2B (Fig. O^c)) 
interpolates between lA and 2n. In the 2B, the lA is split into two sublayers 
consisting of rows that are displaced in height and which can transform smoothly 
into 2n upon increasing the gap. The rhombic phase 27?. (Fig.JSI^d)) is found between 
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2n and 2A. The rhombic phase is also stable between nO and tiA for n < 5, but 
not in the whole region. In addition, we find that at higher n the interpolating 
structures are mainly prism phases. In agreement with experiments, we find two 
types of prisms, one with a square base nVa (Fig. El^e)), and one with a triangular 
base uVa (Fig-Elf)): where n indicates the number of particles in the prism base. As 
shown in Fig. |3Ig),(h), these structures display large gaps as a result of periodically 
repeated stacking faults in the packing which, nevertheless, allow particles to pack 
more efficiently, than a phase consisting of parallel planes of particles. For n > 3, 
differences between simulations and experiments emerge. We find that the stability 
region of interpolating structures between nA, and {n + l)n decreases for larger H, 
becoming invisible on the scale of Fig. [5] for 19 = 5 A -^ GD. On the other hand, the 
region of stability of the interpolating structures between nO and nA increases while 
increasing the wall separation, becoming stable also at low packing fractions for the 
transitions 18 = SD -^ 5A, and possibly 710 — GD — > 6A. We also note that the 
solid-solid transitions are first-order with a clear density jump at low rj, but they get 
weaker (and maybe even continuous) upon approaching the maximum packing limit. 
In addition, the rhombic and buckling phases are highly degenerate as we find zig-zag 
and linear buckling or rhombic phases, and a combination of those. 

We now turn our attention to the fluid-solid transition. In Fig. 01 we plot 
the chemical potential /?//'^^p at the freezing transition of the confined system as a 
function of H . The freezing for crystal slabs with a triangular symmetry are denoted 
by triangles, while the square symmetry is displayed by squares. We find strong 
oscillations in the chemical potential reminiscent to the (in)commensurability of the 
crystal structures with plate separation. The highest values for P^f-^^ are reached 
at the transition region nA ^ {n -\- l)n, corresponding to plate separations where 
both structures are incommensurate and hence unfavorable. In this regime, /3/i'^^P can 
reach values that are higher than the bulk freezing chemical potential /?^''""^ (the black 
vertical line in Fig.Q, corresponding to capillary melting, while the freezing transitions 
with /3/i'^^P lower than the bulk value correspond to capillary freezing. Hence, we find 
a reentrant capillary freezing/melting behavior for wall separations 1 < H /a < 3.5. 
In addition, we compare our results with the predictions of the Kelvin eauation|21| : 
^^cap _ /3^i^"ik _ TT(j^/3H{-fiuf — ^ws)/{ris — Vf) using the parameters determined 
in our simulations. The thick dashed line in Fig. 01 is the prediction of the Kelvin 
equation for the (111) crystal plane (triangular order) at the walls, while the dotted 
line is that for the (100) plane (square order). The Kelvin equation predicts capillary 
freezing for the triangular structure and capillary melting for the square structures. 
The Kelvin equation predictions are in reasonable agreement with our simulations for 
triangular order for wall separation as small as H/a ~ 4, but deviates for smaller H, 
while the prediction for the square structure is in agreement only at very small H. 
It is surprising to find qualitative agreement at small H since the Kelvin equation is 
valid in the limit H/a — > oo. 

4. Conclusion 

In summary, we have calculated the equilibrium phase diagram of confined hard 
spheres using free energy calculations with a novel integration path. The high density 
sequence of structures is in good agreement with experimental results. We find that the 
prism phases are thermodynamically stable also at lower densities, and this work will, 
hopefully, stimulate further experimental investigations, for a quantitative comparison 
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Figure 3. (Color online) Stable solid structures of confined hard spheres, (a) 
The triangular phase 2A (b) The square phase 20 (c) The buckling phase 2B 
(d) The rhombic phase 27?. (e),(g) The prism phase with square symmetry S'Pi:] 
(f),(h) The prism phase with triangular symmetry SPa- Iii (^)"(f) the point of 
view is at an angle of 30° to the z direction. In (g),(h) the point of view is at 
an angle of 90°. Different shades (colors) indicate particles in different planes 
((a)-(d)) or particles belonging to different prism structures ((e)-(h)). 
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Figure 4. Chemical potential /3/i at fluid-solid coexistence, for different waU 
separations Hja. Tfie symbols are the simulation results for the triangular (A) 
and square structures (D). The thin dashed line is a guide to the eye. The 
thick continuous line indicate the value of the bulk freezing chemical potential 
/3/x = 16.08. The thick dashed and dotted curves are the prediction of the Kelvin 
equation for the 111, and 100 planes parallel to the walls, respectively. 



at intermediate packing fractions. In addition, our results show an intriguing sequence 
of melting and freezing transitions upon increasing the distance between the walls of a 
slit which is in contact with a bulk reservoir. The mechanical behavior is therefore very 
sensitive on the degree of confinement, and the knowledge of the phase diagram can 
help the understanding and fabrication of new materials. The transition from confined 
to bulk behavior, and the interface between different solid structures (studied in lower 
dimensions in Ref. |22j ) represent interesting directions for future investigations. 
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